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Abstract 

Different approaches to describe initial stages of relativistic heavy ion collisions 
are discussed qualitatively and quantitatively. An Effective String Rope Model is 
presented for heavy ion collisions at RHIC energies. Our model takes into account 
baryon recoil for both target and projectile, arising from the acceleration of partons 
in an effective field, F^ u , produced in the collision. The typical field strength (string 
tension) for RHIC energies is about 5 — 12 GeV/fm, what allows us to talk about 
"string ropes". The results show that a QGP forms a tilted disk, such that the 
direction of the largest pressure gradient stays in the reaction plane, but deviates 
from both the beam and the usual transverse flow directions. The produced initial 
state can be used as an initial condition for further hydrodynamical calculations. 
Such initial conditions lead to the creation of third flow component. 
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1 Introduction 



Fluid dynamical models are widely used to describe heavy ion collisions. Their 
advantage is that one can vary flexibly the Equation of State (EoS) of the mat- 
ter and test its consequences on the reaction dynamics and the outcome. This 
makes fluid dynamical models a very powerful tool to study possible phase 
transitions in heavy ion collisions - such as the liquid-gas or the Quark-Gluon 
Plasma (QGP) phase transition. For example, the only models that can han- 
dle the supercooled QGP are hydrodynamical models with a corresponding 
EoS. For highest energies achived nowadays at RHIC hydrodynamic calcula- 
tions give a good description of the observed radial and elliptic flows [1-6], in 
contrast to microscopic models, like HIJING [7] and UrQMD [8]. 

In energetic collisions of large heavy ions, especially if a QGP is formed in the 
collision, one-fluid dynamics is a valid and good description for the intermedi- 
ate stages of the reaction. Here, interactions are strong and frequent, so that 
other models, (e.g. transport models, string models, etc., that assume binary 
collisions, with free propagation of constituents between collisions) have lim- 
ited validity. On the other hand, the initial and final, Freeze Out (FO), stages 
of the reaction are outside the domain of applicability of the fluid dynamical 
model. 

Thus, the realistic, and detailed description of an energetic heavy ion reaction 
requires a Multi Module Model, where the different stages of the reaction are 
each described with a suitable theoretical approach. It is important that these 
Modules are coupled to each other correctly: on the interface, which is a three 
dimensional hypersurface in space-time with normal <icr M , all conservation laws 
should be satisfied, and entropy should not decrease. These matching condi- 
tions were worked out and studied for the matching at FO hypersurfaces in 
details in Refs. [9-11]. In Ref [12] it has been shown that the entropy con- 
straint is very sensitive to the pre-FO EoS, i.e. the arbitrary choice of the 
pre-FO state does not necessarily lead to the physical FO. 

Similar ideas form the base of Ref. [13,6,14] where authors combined two 
modules - hydro and UrQMD - replacing the hadronic phase of hydrodynamics 
with hadronic transport model to describe properly chemical and thermal FO. 

The initial stages are the most problematic. We will discuss here different 
models for the initial state in section 2 and will present our modeling of it in 
section 3. 

Our goal is to build a Multi Module Model for ultra-relativistic heavy ion col- 
lisions that is valid for RHIC and LHC energies, and maybe for the most ener- 
getic SPS collisions. The present work is just the first step, but an important 
one: the Effective String Rope model is developed for the most problematic 
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module - module describing the initial stages of collisions. 



2 Initial state of Relativistic Heavy Ion Collision 

The situation at the beginning of an energetic heavy ion collision is highly- 
complicated and not easily understandable. None of the theoretical models 
currently on the physics market can unambiguously describe the initial stages. 

As was pointed out in the introduction, at ultra-relativistic energies the one- 
fluid dynamical models can not be justified for such a situation. Frequently, 
two- or three-fluid models are used to remedy the difficulties and to model the 
process of QGP formation and thermalization [15-17]. Here the problem is 
transferred to the determination of drag-, friction- and transfer- terms among 
the fluid components, and a new problem is introduced with the (unjustified) 
use of an EoS in each component in a nonequilibrated situation where an 
EoS is not defined. Strictly speaking this approach can only be justified for 
mixtures of noninteracting ideal gas components. 

The kinetic transport models (see for example, Ultra-relativistic Quantum 
Molecular Dynamics (UrQMD) [18]) and parton cascades (see for example, the 
Molnar Parton Cascade (MPC) [19,20]) are also frequently applied to describe 
the high energy collisions from the very beginning. Some of the parton cascades 
(for example, MPC) are able to handle only 2^2 reactions so far, what 
cannot be justified for dense and strongly interacting systems. The parton 
cascades also have to approximate media-effects, since the parton structure 
function, measured for a free nucleon, is modified inside the nucleus. This is 
the so-called shadowing effect [21]. This approximation procedure is also not 
unambiguous, see for example Ref. [22,23] for more details. 

Some models, like UrQMD, were initially developed for much lower energies 
than available nowadays, but even then, when the particle propagates in a 
hot and dense and nonequilibrated medium of highly excited hadrons, the 
properties of this particle might change significantly. Thus, properties like ef- 
fective masses, decay widths, effective momenta, and in-medium cross sections 
have to be calculated for a given local situation at the point where particle 
propagates. Unfortunately this is a very complicated and numerically very 
time consuming task, so in most of the models drastic approximations are 
made (see for example [18]). In a relativistic heavy ion collision, especially 
at the early stages, the excitation energy per particle would be so high that 
the concept of resonance (we have to deal not only with nucleons, but with 
many different baryons and mesons) has no sense any more. Here the quark 
and gluon degrees of freedom must be taken into account. Therefore, some 
models, including UrQMD, replace the resonances by continuous string exci- 
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tations for internal excitation energies higher than approximately 2 GeV. The 
properties of such effective "strings" are basically unknown and can be varied 
within a wide range, almost like free parameters [24] . 

Perturbative Quantum Chromo Dynamics (pQCD) would be, in principle, 
the proper model for describing our systems at very high energies. Unfortu- 
nately, pQCD itself is not applicable for heavy ion reactions at RHIC ener- 
gies. Nevertheless, the pQCD calculations with some extra non-perturbative, 
phenomenological assumptions, like saturation of a gluon plasma, can be per- 
formed. Different models following this scenario have been proposed in [25-30] . 
The extension of the pQCD+saturation model for nonzero impact parameters, 
which is of primary interest for us, leads to serious complications in calcula- 
tions [29]. 

Thus, since we can not unambiguously describe the initial stages of relativistic 
heavy ion collisions in microscopical models, we have to use some phenomeno- 
logical models (in fact all the microscopical models discussed above already 
have some phenomenological elements). 

2.1 Different models 

There are two phenomenological models, which are most frequently discussed 
in the literature, and which claim to give us at least a qualitative understand- 
ing of the initial stages of the relativistic heavy ion collisions. 

2.2 The Landau model 

The Landau model has been developed to generalize the fluid dynamical ap- 
proach to the energy range where we can not neglect relativistic effects, such as 
Lorentz contraction of the colliding nuclei. Actually the model was originally 
developed for pp collisions, but the generalization for AA is rather straight- 
forward and the fluid dynamical expansion is even more justified for bigger 
systems. 

It is the most appropriate model for relativistic heavy ion collisions in the en- 
ergy range E[ ab = 10 — 100 A ■ GeV. For energies above r ) CM ~ 3 — 10, so in the 
Center of Mass (CM) frame the colliding nuclei are considerably contracted 
and have disk-like shapes. The Landau fluid dynamical model assumes the 
initial thermalized state as a static, completely stopped, homogeneous disk, 
contracted in the beam direction by the factor 7cm (see Fig. 1). The fluid 
dynamical expansion initially starts in the direction of the largest pressure 
gradient - in the beam direction in our case. Only when the system has ex- 
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Before collision Initial Landau's state 

Fig. 1. Initial state in Landau model. 

panded in the z direction to a size comparable to its transverse diameter the 
transverse expansion also considered. In this second phase the expansion is 
2+1 dimensional, while initially it was 1+1 dimensional (all quantities were z 
and t dependent). 

The inputs to this model are: A) the equation of state; B) the initial energy 
and baryon densities; C) the FO conditions. Usually the initial energy density, 
eo, is estimated in the following way: 



e = E CM /V , where V = V rest /^ 



CM 



2.3 The Bjorken model 



If we increase the collision energy further to Ecu = 100 A ■ GeV or more, the 
expectation is that the concept of complete stopping does not work anymore. 
Furthermore, when two nucleons from the target and the projectile nuclei 
collide at these energies, their momenta are so high that they see each other 
as compound objects, containing three valence quarks and other partons. The 
point is that the proper participants at these energies are the partons and not 
the nucleons. The initial momenta of partons are very high and we expect these 
two parton clouds to cross each other and keep going with almost their original 
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rapidities. Thus, contrary to complete stopping assumed in the Landau model, 
we have now the opposite scenario - nuclei become essentially transparent to 
each other. 



During their interpenetration, however, the target and projectile partons may 
exchange color charges, what leads to the creation of a chromoelectric field in 
the middle of the reaction volume. The classical analogy, that one can keep in 
mind for a better understanding, is the electric field between two condenser 
plates, moving away from each other. From QCD calculations we know that 
such a chromoelectric field should be confined in the transverse direction and 
have the geometry of a string. The energy density of such a string is substantial 
- already for the strings produced in hadronic collisions the string tension 
is about 1 GeV/fm , and the theoretical simulations show that even more 
energetic objects may be produced in heavy ion collisions. Later these strings 
start to decay by producing quarks, antiquarks and gluons, but because of 
baryon number conservation, the matter in the middle of reaction volume is 
baryon free. Theoreticians like very much this last feature of the model, since 
the baryon free matter is rather pleasant for calculations. 



A 



Before collision 




Initial Bjorken's baryonfree "fireball" 



Fig. 2. Initial state in the Bjorken model. The dash-dotted lines denote the chro- 
moelectric field strings. 



The important input to the model is given by experimental data from ultra- 
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relativistic pp collisions. The rapidity distribution of the charged particle mul- 
tiplicity, dN c h/dy, can be divided, with reasonable accuracy, into three re- 
gions. Two smooth peaks around y = ±y , where y^ is the original rapidity, 
are believed to be originated by target and projectile partons, interpenetrat- 
ing through each other almost without rapidity lost. Between those we can 
see the plateau, - dN c h/dy = const, - generated by our baryon free flux tube. 
Assuming the same symmetry, dQ/dy = const, for all other quantities a sim- 
ple hydrodynamical model with an analytical solution has been developed 
by Bjorken in Ref. [31]. The above assumption, dQ/dy = const, means that 
our model is boost invariant - all the quantities are functions of the proper 
time only, Q = Q(r), where r = y/t 2 — z 2 . Obviously the model has a 1+1 
dimensional structure. 

In order to study transverse expansions, the Bjorken model has to be sup- 
plemented with a description of this process separately, see for example [32]. 
The boost invariance of the model is also very much appreciated by theoreti- 
cians, since it is then sufficient to perform simulations for the y = (z = 0) 
plane only, and then boost the result to reproduce all of the phase space. The 
Bjorken model does not tell us anything about the energy density distribution 
in the plane perpendicular to the beam direction, say the z = (y = 0) plane; 
therefore additional assumptions are required. If we define a grid in the [x, y] 
plane, then the energy density profiles proportional to the number of binary 
collisions, c(x i ,y j ,z = 0) oc N target (xi,yj) ■ N target (xi,yj), or proportional to 
number of participants, c(x i ,y j ,z = 0) oc {N target (xi,yj) + N target (xi, %■)), are 
the most popular assumptions in literature (see [33] for recent overview). 

The initial energy density in the Bjorken model can be estimated in the fol- 
lowing way. At the y = (z = 0) plane the full energy is the transverse energy, 
Ef u u = E T . One may also notice that for the Bjorken boost invariant model, 
the pseudo-rapidity (n = — In tan |), where 9 is a polar angle of emitted parti- 
cle, is the same as rapidity, n — y. If in the experiment we measure the energy 
AEt emitted in a very small rapidity window, [0, Ay], then we can estimate 
the initial energy density: 



Now using the following expressions: transverse area A = nR 2 A , where Ra is 
the radius of nucleus with A nucleons, and size in the beam direction Az = 
r An = t A y, where tq is the thermalization time, - we obtain 



Here the thermalization time is unknown and should be estimated from some 





(2) 
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other model for the particular reaction (for example, To 1 fm is usually 
taken for SPS Pb+Pb 158 GeV/nucl collisions). Another question is - can we 
really use the post FO final state measurement, , directly to estimate 

the initial energy density. Briefly speaking, this means to neglect the longitu- 
dinal expansion caused by pressure during the evolution of our fireball, which 
would reduce our initial energy density to the one estimated from the data 
by eq. (2). Therefore, one may expect that the initial energy density is larger 
than the one given by this simple estimate. 



2.4 Bjorken model with baryon recoil 

Certainly, the assumption of non-stopping partons is an idealization: from 
basic physical principles we may expect that if two moving quarks stretch 
a string between themselves with a substantial string tension, they should 
be decelerated, since during such a process they are converting their kinetic 
energy into the energy of the string. Thus, the creation of the chromoelectric 
field should produce some baryon recoil for both colliding nuclei. In the first 
approximation this was done in Ref. [34]. 

The authors developed the model for pA collisions. The initial momentum of 
the projectile, p, was considered to be so huge, that the recoil effect could be 
disregarded for it, while for the target partons the recoil arises from the accel- 
eration of partons in an effective external field F^", produced in the collision. 
This baryon recoil was included in a way that guarantees the conservation 
of baryon current. At the same time the energy- momentum conservation was 
obviously violated - the effective field produced in the collision was considered 
to be external, and the projectile matter was not taken into account in the 
energy-momentum tensor of the system. Thus, the authors ended up with the 
following system of the modified fluid dynamical equations: 

d,T^ = F^N, + , (3) 

d,N» = , (4) 



where T^ v = e t ((1 + Cg) u^u v t — c\g^ v ) is an energy-momentum tensor (only 
target matter is taken into account), c\ comes from the EoS P = c^e, is the 
pion source term [32], and iV M = nw M . The effective field, F^, is parametrized 
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in the following way: 



-a\ 





(TOO / 



(5) 



where a = const is the field strength, or string tension, of the above described 
effective external field. The authors have found o »s 3 — 4 GeV/ fm. 

Without the source term the system (3, 4) can be solved analytically [34]. 
The source term is necessary at the later stages of the reaction, since it allows 
the neutralization process, i.e. the mechanism by which the energy stored in 
the field is converted back into the energy of matter (see [34] for details). 
The interplay between the baryon recoil and the energy, coming from the 
neutralization of the strings - string decay, leads to much higher energy density 
in the fragmentation region than the original Bjorken model predicts. We will 
come up with the modification of this nice simple model in section 3. 

Experiments have now entered the region where the Bjorken model is ex- 
pected to be applicable, but there is no clear and unambiguous confirmation 
that Nature follows this scenario. The preliminary experimental results from 
RHIC do not show transparency - most particle multiplicities do not show a 
dip in the rapidity spectra, but rather a plateau around mid-rapidity[35,36], 
which is a sign of strong stopping. Furthermore, a very strong elliptic flow 
(t>2 flow component^) has been measured, which shows a clear peak around 
mid rapidity [35,40]. To build such a strong elliptic flow, strong stopping and 
momentum equilibration are required. Also the p/p ratio at mid-rapidity mea- 
sured at RHIC [41,42] (preliminary) is still far from one, which tells us that 
the middle region is not baryon-free. 



Nowadays the expansion in Fourier series is usually applied to study the azimuthal 
distribution of particles [37,38]: 

oo 

1 + 2v n cos(n^») , 
n=l 

where (ft is the azimuthal angle. The first term in square brackets represents the 
isotropic radial flow, while the others are referred to anisotropic flow. The first 
Fourier coefficient v\ = (cos (ft) = (p x /pt) is called directed flow, and the second one 
V2 = (cos (2(f))} = {(px/pt) 2 — (Py/Pt) 2 ) is called elliptic flow. 
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2.5 Third flow component 



As discussed in the introduction, the fluid dynamics is governed by the Equa- 
tion of State of the matter, and the analysis of the resulting flow patterns 
turned out to be one of the best tools to extract the EoS from the outcome of 
a heavy ion collision. The final event shape must carry the information about 
the pressure development during the collision including the early stages of the 
collision. 

The phase transition to the QGP is connected to a decrease of pressure ac- 
cording to most theoretical estimates, not only in strong first order phase 
transition models, but even if we have a smooth but rapid, gradual transition. 
This reduced pressure around the phase transition threshold are known for 
long (see e.g. [45]), and it was emphasized as a possible QGP signal: this "soft 
point" of the EoS might be possible to observe in excitation functions of col- 
lective flow data [48,49]. It was pointed out earlier [43] that the decrease in the 
out-of-plane (squeeze out) emission is even more sensitive to the pressure drop 
and it decreases due to plasma formation stronger than the in-plane collective 
flow. 

In this section we want to discuss another consequence of the same softening 
in the EoS, which is a recently (1999) identified new distinct flow pattern 
- the third flow component [44]. This flow pattern can be seen in almost 
all theoretical fluid dynamics calculations with QGP formation, and also in 
experimental data, but it was not discussed earlier, see for example [43,46- 
48,50]. 

Two basic flow patterns have been predicted and detected before at lower 
energies: 

A) the directed transverse flow in the reaction plane, or side-splash, or bounce 
off, which is most frequently presented on the well known p x vs. y diagram (or 
Vi flow component vs. y) and seen at all energies in heavy ion collisions from 
energies of 30 A ■ MeV to 165 A ■ GeV [54-57]; 

B) the squeeze-out effect which is an enhanced emission of particles transverse 
to the reaction plane at center-of-mass rapidities. 

At lower energies the directed transverse flow resulted in a smooth, linear p x 
vs. y dependence at CM rapidities. This straight line behavior connecting the 
maximum at y pro j and the minimum at y targ was so typical that it was used 
to compare flow data at different beam energies and masses. 

At higher energies, say above 10 — 11 A- GeV in Lab frame, the deviations from 
this straight line behavior have been observed [51,54,56,57]. Unfortunately, 
directed flow has not yet been measured at RHIC. 
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In theoretical calculations the important feature of this third flow component 
is that it is clearly seen in fluid dynamical calculations when, and only when, 
the QGP formation was allowed during the calculation (see [44] for overview). 
In sharp contrast, the solutions with a purely hadronic EoS did not show this 
effect, and the maximum and minimum of the p x curve could be connected 
with a straight line[f]. 



2.6 Possible source of the third flow component in fluid dynamics 



V 



Before collision 




Tilted initial Landau's disk 



Fig. 3. Tilted initial state. The direction of the largest pressure gradient stays in the 
reaction plane, but deviates from both the beam and the usual transverse flow direc- 
tions. Such initial conditions may lead to the creation of the third flow component 
[44]. 

The following possible source of the third flow component has been proposed 
in Ref . [44] . If QGP is formed in the collision, strong and rapid stopping may 
take place, and a behavior close to that of a one-fluid model one might be 
established. I.e., let's assume that Bjorken scenario is not realized, at least in 



2 Note that all these fluid dynamical calculations were done much before the ex- 
periments. The first quantitative flow predictions [46] preceded the experiments by 
six years and gave rather good agreement with the data. 
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its form discussed in previous section, but that the Landau fluid dynamical 
model is applicable for collisions of massive ultra-relativistic heavy ions (see 
the discussion at the and of section 2.4). The soft and compressible QGP forms 
a rather flat disk which is at rest in the CM system (see Fig. 3). Then this 
disk starts to expand in the direction of the largest pressure gradient. At small 
but finite impact parameters we may assume that this disk is tilted, and the 
direction of fastest expansion will stay in the reaction plane, but will deviate 
from the both beam axis and the usual transverse flow direction. Thus, a third 
flow component might develop from the large pressure gradient and tilted and 
strongly Lorentz contracted initial state. Now the flow in the reaction plane 
will have an "elliptic structure" [56,58,59] and must be characterized by two 
flow components: the usual transverse flow and the one orthogonal to it - the 
third flow component or antiflow. At the same time as the primary directed 
flow is weakened by the stronger Lorentz contraction at higher energies, this 
third flow component is strengthened by increased pressure gradient arising 
from the Lorentz contraction. 

The concept of a tilted "fireball" was proposed initially in the so-called "fire- 
streak" models [60,61]. The idea behind it is rather simple and has to do with 
momentum conservation. Without going into mathematical details we can un- 
derstand it in the following way. First, let us create a grid in the [x, y]-plane 
and subdivide the projectile and target into streaks parallel to the beam di- 
rection. We try to describe the collision with finite impact parameter. In the 
CM frame in one streak the pieces of matter from target and projectile will 
have the same absolute values of rapidities, but different momenta, since they 
in general contain different amount of matter (see Fig. 4). If the collision hap- 
pens so fast that there is no interaction between streaks, then the momentum 
conservation does allow complete stopping, i.e., the final streak - the piece of 
matter within one tube with the same transverse coordinates after collision of 
target and projectile fragments - will have a momentum defined by momentum 
conservation, and correspondingly its center of mass has to move with some 
rapidity. Crossing the participant region along x axes in the reaction plane 
we will see that in the middle the final streaks rapidities are small, since the 
target and projectile fragments were almost equal; while at the edges of the 
participant region the final streaks move rather fast (see Fig. 4). Thus, after 
some time our tilted "fireball" will degenerate into a "firestreak" . Please no- 
tice the difference between "firestreak" and the tilted Landau disk discussed 
above: the tilted disk was at rest in CM frame, i.e., the complete stopping 
has been assumed, while "firestreak" geometry is dynamic. By "dynamic" we 
mean that different parts of "firestreak" move with different rapidities. 

In "firestreak" models such a geometry has been assumed at the freeze-out 
time, and then, based on energy, charge and baryon number conservations (we 
already conserve momentum) the authors tried to reproduce the parameters 
of the post FO distribution [61]. This procedure is somewhat similar to the 
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Fig. 4. The "firestreak" initial state. 

statistical production model, which is used nowadays to calculate production 
of different hadronic species at chemical FO hypersurface. Such calculations 
were not really supported by experiment. Nevertheless, in our opinion, this 
scenario might be a good qualitative description for the initial stages of the 
heavy ion collisions until the local thermalization achieved. When the matter 
is locally thermalized and the pressure is built up, we cannot neglect hydro- 
dynamical expansion, which will smear out this initial distribution producing 
a more or less spherical fireball before the freeze out process starts. Neverthe- 
less, such a hydrodynamical expansion will initially start in the direction of 
the largest pressure gradient (see Fig. 4) and, together with the initial velocity 
field, might produce a two component (elliptic) structure of the flow in the 
reaction plane, identified in ref. [44]. 



3 Effective String Rope model 

In this section we will discuss so called Effective String Rope model (ESRM) 
for energy, pressure and flow velocity distributions at the beginning of ultra- 
relativistic heavy ion collisions. The output of this model can be used as an 
initial condition for further hydrodynamic calculations. 

One important conclusion of heavy ion research in the last decade is that 
standard 'hadronic' string models fail to describe heavy ion experiments. All 
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string models had to introduce new, energetic objects: string ropes [62,63], 
quark clusters [64], or fused strings [65], in order to describe the abundant 
formation of massive particles like strange antibaryons. Based on this, we de- 
scribe the initial moments of the reaction in the framework of classical, or 
coherent, Yang-Mills theory, following the Bjorken model with baryon recoil 
Ref. [34] (see section 2.4) assuming a larger field strength (string tension) 
than in ordinary hadron-hadron collisions. For example, calculations both in 
the Quark Gluon String Model (QGSM) [66,46,67] and in the Monte Carlo 
string fusion model [65] indicate that the energy density of strings reaches 
8 — 10 GeV/ fm already in SPS reactions. This is nearly 10 times more than 
the tension used in standard, 'hadronic', string models where a ~ 1 GeV j fm. 
In addition we now satisfy all conservation laws exactly, while in Ref. [34], as 
it was discussed in section 2.4, energy-momentum conservation was irrelevant 
and was ignored. Thus, in this approach for the first time the initial trans- 
parency/stopping and energy deposited into strings and "string ropes" was 
determined consistently with each other [70] . 

We do not describe the hadronization of our strings. Thus, this approach can 
be justified only for the very initial stages. The effective way to perform string 
decay will be discussed in section 3.3. 

Recent parton kinetic models indicate that saturation of a gluon plasma takes 
place in a very short time - r sat = 0.09,0.27 fm/c for LHC, SPS energies 
respectively [27], and a pressure needs r p = 1 — 2 fm/c (SPS) and t p = 3 — 
5 fm/c (LHC) to be established [39]. More importantly the first experiments 
at RHIC yield strong elliptic flow, which, as it was mention in introduction, 
was not reproduced in any other model, except in fluid dynamical models 
with a QGP EoS [35]. This is a strong experimental indication that transverse 
pressure builds up early in these reactions, in a few fm/c, and strong stopping 
is also necessary to create strong flow before Freeze Out, which usually happens 
when the radius of the system is not more than 10 fm. Our model builds up the 
locally equilibrated matter - the initial conditions for further hydrodynamical 
expansion - at ti a b = 3/5 fm / c for central / peripheral regions at RHIC energies, 
which is in agreement with previous estimations as well as the data. 

The presentation in this section is based on Refs. [68-70], but we want to stress 
that the preliminary results were presented at conferences and published in 
their proceedings, while the model was still being developed and, thus, Ref. 
[70] is more up to date than Refs. [68,69]. And in this work we also improved 
the presentation of our model, particularly in sections 3.3, 3.4, 3.5. For the 
continuous description we will repeat in this section some details from Ref. 
[70], but we will also underline the differences. 
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3. 1 Formulation of model 



The basic idea is to generalize the Bjorken model with baryon recoil, presented 
in section 2.4, for collisions of two heavy ions and improve it by strictly satis- 
fying conservation laws [68-70]. Similar ideas were also presented recently in 
[71]. 

First of all, we would like to create a grid in the [x, y]-plane (z is the beam 
axes, [z,x] is the reaction plane). We will describe the nucleus-nucleus col- 
lision as a sum of independent streak-by-streak collisions corresponding to 
the same transverse coordinates, {xi,yj}. We assume that baryon recoil for 
both target and projectile arise from the acceleration of partons in an effec- 
tive field, F^, produced in the interaction. Of course, the physical picture 
behind this model is based on chromoelectric flux tube or string models, but 
for our purpose we consider F^ u as an effective Abelian field. The most im- 
portant consequence of the non- Abelian fields, i.e., its self interaction and the 
resulting flux tubes of constant cross section, are, nevertheless, reflected in our 
model: we assume that the field is one-dimensional. The fields generated by 
the colliding streaks are of constant cross section during the whole evolution, 
and only their lengths increase with time. As the string tension is constant, 
the energy of the string increases linearly with increasing length. The single 
phenomenological parameter we use to describe this field must be fixed from 
comparison with experimental data. 

We describe the streak-streak collision based on the following set of equations: 

^E^ = E^V\^, (6) 

i i 

c> M ]T7Vf = 0, ? = 1,2, (7) 

i 



where iVf is the baryon current of ith nucleus, is the color charge (discussed 
later in more detail). We are working in the Center of Rapidity Frame (CRF), 
which is the same for all streaks. We will use the parameterization: 

= riiU^ , uf = (cosh?/;, sinhyj) . (8) 



J 7 *" is the energy-momentum tensor. It consists of five parts, corresponding 
to two nuclei and free field energy (also divided into two parts) and one term 
defining the QCD perturbative vacuum. (While in the similar equation in 
the Bjorken model with baryon recoil, eq. (3), contains only the target 
matter.) 

/ j ^ i ' -*■ pert 
i=l,2 
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E h (( X + c 9 « - c o^) + ^] + BgT ■ (9) 

i=l,2 



Here B is the bag constant, the equation of state is Pi = c^ei, where and Pi 
are the energy density and pressure of QGP. 

Comparing our equation (6) with the corresponding one in the Bjorken model 
with baryon recoil, eq. (3), one may notice that we neglected the pion source 
term. As it was discussed in section 2.4 the source term takes care about 
the string hadronization process. We are planning to use our model at the 
beginning of the ultra-relativistic heavy ion collisions only, and the string decay 
process at later stages will be described separately. Section 3.3 is devoted to 
this question. The advantage of equations (6, 7) is a simple analytic solution. 

Within each streak we form only one flux tube with uniform field strength 
or field tension, ex, between the target and the projectile streaks. The field 
extends to the outside edges of the target and projectile streaks, i.e., the field 
also overlaps with the space-time domains of the matter streaks. For practical 
and symmetry purposes we, however, divide this field into two spatial domains, 
a target and a projectile domain, (i = 1,2), separated at a fixed point, z sep , 
so that (Ji = <r 2 = <y ■ The choice of this point will be specified later. (The field 
is constant and the only change is that it extends with time at its two ends.) 

In complete analogy to electro-magnetic field 

/ 



F? v = d^A" - dTAf = 



-Oi 
(Ji 



(10) 

\°i o J 

a, = <9 3 A° - d ^ , (11) 



dCp ■ 
p <9(<9Mf) 

£f,i = —^Fi^Ff . (13) 



In our case the string tensions, <7j, will be constant in the space-time region 
after string creation and before string decay. The creation of fields will be 
discussed later in more detail. 

Since during the interpenetration both nuclei are closely following the light 
cone, the geometry of the reaction suggests to use the light cone variables 

(z,t) -> (x + ,x~), x ± = t±z, (14) 



16 



and to assume that ei,yi,ni are functions of x~ only (the first nucleus goes 
from right to left, yi(0) = — y < 0) and e 2 ,y2,n 2 depend on x + only (the 
second nucleus goes from left to right, 1/2(0) = yo > 0). 

In terms of lightcone variables (for the transformation matrixes between space- 
time and lightcone coordinates see Ref. [70]): 



, f - //,(//'' : <r') = 11 ,( : ?y , 

h i+ e 2y * hi_ 
y hi_ h i+ e~ 2yi 

— c )ej . 







\T- + T--j 


n 


where 




hi+ = (1 + cDd , 


hi- = 



FA 



(15) 
(16) 



(17) 



The other tensors in light cone variables are: 



Fi = 



-pert 



(1 


7++ 

i 


Ft 


w 


F~ 






:) 









Oi 



-<Ji 



(18) 



(19) 



The energy-momentum tensor for free field in the light cone variables is: 



1 K 2 

a? 



(20) 



At the time of the first contact of the two streaks, t — 0, there is no string 
tension. We assume that strings are created, i.e., the string tension achieves 
the value a at time t = to, corresponding to complete penetration of streaks 
through each other (see Fig. 5). 



3.2 Conservation laws and string creation 



In lightcone variables eq. (7) may be rewritten as 
d-Ni + d+N+ = . 



(21) 
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Fig. 5. Streak-streak collision. The time of first touch of streaks is t = 0, and 
t = to corresponds to complete penetration of streaks through each other. At 
this time strings are created, i.e., the string tension reaches an absolute value 
a = A (^) 2 no^/lih (32). The dash-dotted lines denote the chromoelectric field 
strings (it will be shown later that a is so big that we can talk about several paral- 
lel strings or a "string rope"). 

So, we have a sum of two terms, each depending on different independent 
variables, and the solution can be found in the following way: 

<9_iVf = a, d + N 2 + = -a , 

1 + 2 (22) 

jVf = ax~ + (JVf ) , N+ = -ax + + (N+) , 

where the index indicates the initial proper density, which is the normal 
nuclear density: no = 0.145 fm~ 3 . Since both iVf and iV^ are positive (and 
also more or less symmetric) we can conclude that for our case a = 0. 

Finally, 

JVf = n ie ~ yi = n e yo , iV 2 + = n 2 e V2 = n e yo , (23) 
n x = n e yo+yi , n 2 = n e y °- y2 . (24) 



where y (—yo) is the initial rapidity of nucleus 2 (1) in the CRF, respectively. 
The other components are given by eq. (15). 

Let us make an analogy to the electro-magnetic field, where two charges, q x and 
— q2, move in the opposite directions, creating string-like field between them, 
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E = (0, 0, E), which is constrained transversally into a constant cross section. 
The Z-axis goes through the charges, q± and —q 2 , and directed from q± to — q2 
(let us assume that we have such a field configuration). So, forces acting on 
our charges, q\E and —q 2 E, have opposite sign and both are working against 
the expansion of the "string" . In our effective model we use color charges, and 
assume that the vectors of these color charges point in the opposite directions 
in the color space [65], so that the forces acting on both target and projectile 
partons are opposite, both stopping the expansion of the streak. As our field 
strength (string tension, a) is not yet defined we normalize the charges to 
unity: 

qi = —q 2 = 1 > , while = cr 2 = cr . (25) 



Then, we have the forces acting in z direction: q\0\ = a , and q 2 a 2 = —cr. 
Notice again that after string creation fields cr\(x) and a 2 (x) are spatially sepa- 
rated as well as the baryon densities, ri\ and n 2 ; i.e., after complete penetration 
of the initial streaks through each other (see Fig. 5), cr 2 acts on the partons 
on the right side of the separating point z sep = (Z x — l 2 )/2, while the <j\ acts on 
those on the left side. In the absence of matter, in the middle both fields are 
identical, so the exact position of the separating point does not play any role 
until it does not enter the target or projectile matter. The fields <7j are gener- 
ated by the corresponding four-potentials A it which are different and spatially 
separated in the same way (see the dotted separating line, z sep = {l\ — Z 2 ) /2, 
in Fig. 5). 

As described above we do not generate the chromoelectric field self-consistently 
as a product of color currents, which are also affected by the field. Our effective 
fields are external with respect to the colliding partons; that's why we can 
use the expression (20) for the field energy. On the other hand, if we want 
to satisfy the conservation laws, we must generate our effective fields in the 
collision transferring energy from matter to field. It's possible to define new 
conserved quantities based on eq. (6). Using the definition of F^ v , eq. (11), we 
can rewrite eq. (6) as 

d,T^ = £ FTQiN^ = J> [&> (A"N itfi ) - 

i i 

-AV&N^ - d v {AtN^) + Atd»N^\ . (26) 



The solutions for iVf and iV^ 1 ", eq. (23), show that the second term vanishes. 
The fourth term is a vector (Aid^Nf, —A 2 d+N 2 ~) in lightcone coordinates. 
So, if we impose the conditions 

A- = 0, A+ = (27) 
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we can define a new energy-momentum tensor T^ u , such that 

dftv = , (28) 

nniw — \ ^ nh I 7^ _ 

J — 2-^i » pert — 

E (^r - ftA- jvt + ir^Nia) + b<t (29) 



One may notice that the above defined new energy-momentum tensor is not 
symmetric. We know that the energy-momentum tensor is not defined in 
unique way - if T^ u satisfies the eq. d^T^ = 0, then any other tensor T^ v + 
d K ^/^ K , satisfies it also, if = —^/ U ^ K [72]. How to make it symmetric will 
be discussed later in this section. 

To satisfy the above choice of the fields, (25), and imposed conditions, (27), 
we take the vector potentials in the following form: 

A7 =0, At = -a ± x + = -ax + , 

(30) 

A7 = cr 2 x~ = ax~ , A\ = . 



In principle the vector potentials are not defined in unique way - the general 
expression would be: 

A7 = 0, At = —ax + + aire - , 

1 (31) 
A 2 = ox + a 2 x + , A 2 = . 



In our published work, Ref. [70], the choice a\ = a 2 = has been used. It 
appears that this intuitive choice has a deep meaning - with any other a\ and 
a 2 the symmetrization procedure, which will be discussed below, can not be 
applied. 

How to define the a itself? In our calculations we used the parameterization 
(see Refs. [68-70]): 

a = A(^)\^ = Am^, (32) 



where m is the nucleon mass, Eq is the initial energy per nucleon, li and 
l 2 are the initial streak lengths in the CRF (see Fig. 5), and n 2 are the 
numbers of baryons in the initial streaks, Ax A y is the cross section of the 
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streaks, and we've used the relation 7 = eo/m. Similar (ni, n 2 )-dependence 
of a = has been chosen in Ref. [71]. We are working in the system where 
H — c — 1, so a has a dimension of length -2 = energy /length. The typical 
values of dimensionless parameter A are around 0.06 — 0.08 (these values of 
A are for a measured already in GeV/fm). Notice, that there is only one free 
parameter in parameterization (32). The typical values of a are 4—10 GeV/ fm 
for Eq = 65 GeV per nucleon, and a ~ 6 — 15 GeV/ fm for Eq = 100 GeV 
per nucleon. These values are consistent with the energy density of the non- 
hadronized strings in a reaction volume, or "latent energy density", which is 
about 15 GeV/fm 3 for ^ = 200 GeV/nucl [46,66,67]. Our a given by eq. 
(32) grows linearly with Eq. This may be a too strong assumption, since we 
predict much bigger field strength for LHC then the Monte Carlo string fusion 
model [65]. It is easy to modify expression (32) by changing the power 2 to 
1 + a, where a < 1, and probably this will be done after the comparison with 
experimental data. In Ref. [71] a = 0.3 has been chosen based on the low-x 
behavior of the nuclear structure function or parton density [73]. Since such 
an improvement is straightforward, we will keep it in mind, but for further 
discussion we will use the definition (32) as in Ref. [68-70]. 

Using the exact definition of Af, eqs. (30), eqs. (16, 19, 20, 25, 29) and trans- 
formation matrixes from Ref. [70] we obtain 





1 

+ 2 



( h 2+ e^ h 2 _ 
h 2 - h 2+ e- 2 ^ 





' a 2 2b\ 


M 






K 2B a 2 ) 



(33) 




ax~N 2 + 
-ax~N 2 



Notice that the perturbative vacuum term, B, and the free field energy, 
cover all the interacting volume, while the energy densities of matter and 
baryon currents are separated in space. 

Now we are coming back to the question: how to symmetrize T^ u ? One may 
notice that the tensor 



n 



ax'Nl 
ax + N[- 



(34) 
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satisfies the equation d^T^ v = (it's easy to check using the eqs. (23)). Thus 
we can redefine T^ v in a symmetric way: 



rpfll/ rp/XU rpj- 



(35) 



This important step has not been done in Refs. [68-70] . It causes some changes 
in the calculations of the energy densities at t — to (see A), but we will see 
that this does not affect appreciably the final results. 



Now the new conserved quantities are 
Q° = J T 00 dV = Ax A y J f 00 dz , 
qs = J f03 dy = AxAy J f03 dz _ 



(36) 
(37) 



We can rewrite the energy-momentum tensor in (t, z) coordinates: 

+ B 



( J>00 J>03 \ / cr 2 



^30 ^33 



+ 



V 







B 



( 



(ei + Pi) cosh 2 yi — Pi (ei + Pi) cosh y-y sinh y\ 



^ (ei + Pi) cosh yi sinh y x (a + Pi) sinh 2 y x + Pi 



(e 2 + P 2 ) cosh 2 y 2 - P2 (e 2 + P2) cosh y 2 sinh y 2 



^ (e 2 + P 2 ) cosh 7/2 sinh y x (e 2 + P 2 ) sinh 2 y 2 + P 2 y 



+ 



(38) 



ax 



-N 2 ~ N 2 ) 



Obviously, in the absence of the fields, before string creation and after string 
decay, f» v = T» v and, thus, (Q°, Q 3 ) reverts to (P°, P 3 ) - components of the 
four-momenta of the system. 

Based on conservation of Q°, Q 3 we can calculate energy densities, ei(to), e 2 (to), 
at the moment t — t , when the string with tension a is created - eqs. (A. 13, 
A. 14). Another result from A we will need for later discussion is the rapidity 
of the CM frame: 



tanh y CM = Q 3 /Q° = v /M , 



(39) 
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where M = ^±fi . 

Thus, for x ± > x we should solve eqs. (28), with boundary conditions 

Nf(x~ = Xq) = n e Ty " N±(x+ = x ) = n e ±yo 

h 1+ (x~ = x ) = ei(* )(l + cl) h 2+ {x + = x ) = e 2 (t )(l + eg) 

Vi(x~ = x ) = -y y 2 (x + = x ) = y (40) 

a±(x~ = Xo) = a o- 2 (x + = x ) = a 

qi {x~ = x ) = 1 q 2 (x + = x ) = -1 , 



where x = 2t — \z(0)\ defines the string creation surface t = t , for a parton 
or cell element in the position z = z(0) at the time t = 0. The complete 
analytical solution can be found in Ref. [70]. 



3. 3 Hadronization of the String Ropes 



If we let the partons (or cell domains) evolve according to trajectories derived 
in our model (see eqs. (37, 38) in Ref. [70]), the partons from the colliding 
initial streaks will keep going in the initial direction, gradually slowing down 
up to some time t = t^ tU rn] then they will turn and go backwards until the two 
streaks again penetrate through each other (see Fig. 6). Of course, it is hard 
to believe that such a contraction of the reaction volume could really happen 
in heavy ion collisions, because of string decays, string-string interaction, in- 
teraction between streaks and other reasons, which would be very difficult to 
take into account. As we already pointed out earlier we did not install any 
string hadronization mechanism in our model and, thus, there is no way to 
decrease the energy supplied in the strings. If we try to describe the evolution 
of the system with our model after one or both initial streaks turn back, we 
run into the trouble, since the length of the string starts to decrease, so the 
field energy should also decrease, but it has no chance to do this. Therefore 
the above simple analytical model fails to describe the Yo-Yo motion for the 
initial streaks - it has the realistic character until the length of the final streak 
has started to decrease, i.e., approximately up to the dash-dotted line in Fig. 
6. 

In Refs. [68-70] we assumed that the final result of the collision of two streaks, 
after stopping the string's expansion and after its decay, is one streak of the 
length Alf with homogeneous energy density distribution, e/, moving like 
one object with rapidity yj. We assumed that this is due to string-string 
interactions and string decays, which we are not going to describe in our simple 
model. As mentioned above the typical values of the string tension, a, are of 
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Fig. 6. The typical trajectory of the initial streaks with lengths 1% (coming from the 
right) and I2 (coming from the left). 

the order of 10 GeV/ fm, and these may be treated as several parallel strings. 
The string-string interaction will produce a kind of "string rope" between our 
two streaks, which is responsible for final energy density distribution. 

One of the simplest ways to quantitatively take into account string decays is 
presented in Ref. [71]. One of the complications appearing in this case is that 
the hadronization of the strings has to be described in the proper time (i.e. 
all quantities became functions of [t,i]), t = \Jt 2 — z 2 ) which is, in general, 
different for all the pairs of colliding streaks. Then in order to be able to 
transform all of the final streaks to one frame, say CRF, we would have to 
know the details of the their dynamics for a long period of the proper time, 
t = \J t 2 + z 2 . Unfortunately this is not the case in real modeling, for example 
in Ref. [71] authors do not describe the dynamics after all the chromoelectric 
strings for the particular streak pair decay. 

The real situation may be more complicated: when the energy accumulated 
in the strong color fields will be finally released in a production of qq pairs 
and gluons, this process may noticeably change the composition of matter as 
compared to the chemical equilibrium case [74] . Therefore, the matter created 
after the mutual stopping of interpenetrating streaks can not, in general, be 
described by the equilibrium EoS. 

For the stopped and equilibrated matter a set of homogeneous final streaks 
is the simplest assumption. Its advantages are: simple expressions for ej, yj 
and a very simple way to generate an initial state for further hydrodynamical 
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description. We have to give only four numbers (coordinates of the ends of 
the final streak, Z[ e f t and z right , ej and yf) for each point of a transverse grid, 
{xi, Dj} . Later in section 3.5 we present an analytical solution for the expand- 
ing final streaks which seems to be more realistic. Remember - we describe the 
initial state which is not directly observable in experiments. Thus, even a flat 
initial rapidity distribution may end up after the hydrodynamical evolution in 
both a forward-backward peaked or a centrally peaked distribution depending 
on several circumstances. 
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Fig. 7. A) The typical trajectory of the ends of two initial streaks corresponding to 
the numbers of nucleons n\ and 712, £0 = 65 GeV/nucl, A = 0.09 (the parameter 
A was introduced in eq. (32)). Stars denote the points where = yf. From t = to 
until these stars, the streak ends move according trajectories derived in Ref. [70] 
(see eqs. (37, 38)). Then the final streak starts to move like a single object with 
rapidity yj, eq. (41), in CRF. 

B) The same situation as in subplot (A), but for expanding final streak. C) Shows 
e(z) profiles of expanding final streak (t^ = 5 fm). We can clearly see three regions 
- two of forward and backward rarefaction waves and a middle where the initial 
energy density, (eq. (42)), is still preserved. D) Shows the rapidity profile of 
the expanding final streak (th = 5 fm). We can clearly see three regions - two of 
forward and backward rarefaction waves and a middle where the initial rapidity, yj 
(eq. (41)), is still preserved. 

The final energy density and rapidity, ej and yf, should be determined from 
conservation laws. The conservation of the energy and momentum gives for 
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the final rapidity: 



(M 2 (l + c 2 )-2c 2 t, 2 ) + 
cosh 2 y f = ^(M 2 (l + c§) - 2c^ 2 ) 2 + 4c> 2 (M 2 - v 2 ) , (41) 

2(l + c 2 )(M 2 -^ 2 ) 



where we neglected BAlf next to Q / Ax Ay and use notation M, introduced 
in eq. (A. 8); v = tanh|/ is the initial velocity^} Then the expressions for e/ 
is: 



f ((l + 4)cosh 2 y f -cl)Al f ■ 1 J 



The typical trajectory of the streak ends is presented in Fig. 7 (A). From 
t = t they move according to trajectories derived in Ref. [70] (see eqs. (37, 
38)) until they reach the rapidity yi = yf. Later, the final streak starts to 
move like one object with uniform rapidity, yf, until we reach the time when 
the fluid dynamical calculation starts. The time and position of final streak 
formation can be found from the condition y% = yf (see eq. (43) in Ref. [70]). 
We shall call this scenario the "nonexpanding final streak assumption". 

Unfortunately this simple model of the final state causes a principle problem. 
Comparing expressions (41) and (39) one may notice a strange thing: the final 
streak moves like one object with a certain rapidity calculated from energy- 
momentum conservation, but this rapidity, yf, may differ, in general, from the 
rapidity of the Center of Mass of the system given by eq. (39). Obviously, if 
all the elements of the system move with the same rapidity, then the Center of 
Mass of the system also moves with the same rapidity. This problem have been 
pointed out in the Refs. [68-70], but the above oversimplified assumptions have 
been kept because of the simple analytical solution for the e/, yf. We think 
that now we have a deeper understanding of the source of the problem and 
we would like to update the discussion in Ref. [68-70]. This is illustrated and 
elucidated in next section. 



3.4 Flow velocity field generated by pressure in relativistic fluid mechanics 



Let us propose to the reader the following problem: 



3 Please notice the misprint in equation for cosh 2 yf in Refs. [68-70]: (M 2 (l + Cq) + 
2clvl) instead of (M 2 (l + eg) - 2c 2 1 v 2 ). 
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Fig. 8. The homogeneous peace of matter with energy density e and EoS P = P{e) 
moving with velocity vq in the Lab frame. 

Imagine the homogeneous piece of matter with energy density e and Equation 
of State P = P(e) moving with velocity vq in the Lab frame (see Fig. 8). The 
question is to find the energy-momentum tensor, T^ v , and four-momentum, 
P M ; in the Center of Mass frame and in the Lab frame and show that they are 
consistent with each other. 



3.4.0.1 CM frame: In the CM frame, which is in our case the rest frame 
of matter the energy momentum tensor and four-momentum of the system 
are: 



rpfMV 



e ^ 

P 

P 

00 OP; 



(43) 



P M = (E, 0, 0, 0) 



(44) 



where E = eVo, Vq is the volume of our matter in the CM frame, /i, v 



0,1,2,3. 



The four-momentum is connected to the energy-momentum tensor in the fol- 
lowing way [72]: 

P" = [ T^ v dS v , (45) 



where the integral is over any 3-dimensional hypersurface covering all the 
volume, i.e., intersecting all the world lines of the matter in space-time. If we 
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integrate over the hyperplane x° — t — const, then 
P» = J T^dV . 



(46) 



It's easy to check that equations (43) and (44) are indeed connected by the 
relation (46) for any EoS. 



3.4.0.2 Lab frame: In the Lab frame our system is moving with velocity 
v = (0, 0, i>o). The four- momentum undergoes the Lorentz transformation: 



P£ = (P 7o , 0, 0, P 7 o^o) 



(47) 



The invariant expression for the energy-momentum tensor in any reference 
frame is 



= (e + P)u»u u - Pg 



(48) 



In our case = 7o(l, 0, 0, v ). Thus, 

/ ( e + P) 7o 2 -P (e + P^vo \ 



rpflV 



P 

OP 

V (e + P) 7o S> ( e + P)7X + Py 



(49) 



Now we have to use eq. (46) to find P£ from T£ v . 



dV 



e + P)7 2 -P = (e + P) 7o 2 -P V / 7 o , 



Pi = J dV(e + P)>y*v = (e + P)-yfroH/7o , 



(50) 
(51) 



where Vo/70 is the volume of our system in Lab frame. Comparing eqs. (50, 
51) with eq. (47) we obtain the following system of equations: 



(e+P) 7 g-P 



70 



= e7o 



(e + P)7oW = e'jovo 



(52) 



and we see that for a consistent description we have to put P = 0. Why do 
we get such a strong restriction on our EoS?! 
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Fig. 9. All the matter on the hypersurface t = to in the Lab frame has the same 
velocity, but in the CM frame on the hypersurface t± = t® different parts of matter 
have different velocities. 

Certainly the problem is that we made wrong assumptions. We assumed that 
all the parts of our system are moving with the same velocity. This is true 
for the dust, where the pressure is exactly zero. If P > such a situation 
may happen for some particular t in the particular reference frame (remember 
we integrate T^ u over the hypersurface t = const). If we boost our system 
the parts of matter on the new integration hypersurface t\ = const will not 
move with the same velocity, as illustrated in Fig. 9: in the Lab frame all 
the matter on the hypersurface t = to has the same velocity, but in the CM 
frame we have to either integrate over hypersurface t% = t®, where different 
parts of matter have different velocities or integrate over hypersurface t = to 
homogeneous in velocity, but in this case dS^ ^ (dV, 0) in eq. (45). Briefly 
speaking the presence of pressure requires, in general, the flow velocity field 
in CM frame, VcM{t,r), for real physical situation, and our invariant energy- 
momentum tensor, T^, denned by eq. (48) seems to be clever enough to know 
that. 

Thus, as we see, the assumption of the final streak moving like one object 
oversimplifies the real physical situation. The way to improve the model is to 
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let the final streak expand with the rapidity field which will allow us to satisfy 
exactly all the conservation laws. We can describe such a streak expansion 
based on the analytical solution for one-dimensional expansion of matter into 
the vacuum. 



3. 5 Expanding final streaks 



Apart of the principal reasons discussed above, our test results [79,80] indi- 
cated that the uniform final streaks without expansion lead to somewhat un- 
realistic expansion patterns. Thus we had to improve the model and include 
the expansion of the final streaks. 

In this section we discuss the one-dimensional expansion of the finite streak 
into the vacuum (generalizing the description in [81]). The initial condition is 



e(z,t ) 



v(z,t ) 



- oo < z < R- , 
e , R- < z < R + 

, R + < z < oo , 

— 1 , — oo < z < R- 
tanhyo ) R- < z < R + 

1 , R + < z < oo , 



(53) 



(54) 



where R + are the borders of the system. The choice v = ±1 in the vacuum 
is purely conventional, but it guarantees a continuous hydrodynamic solution 
at the boundary to the vacuum, since in the limit of infinite dilution the 
velocity of matter approaches unity. 

Now consider a simple wave, i.e., continuous hydrodynamic flow. This flow re- 
spects entropy conservation [82]. For isentropic flow, the t— and z— component 
of the equation of motion, d^T^ = 0, can be combined to yield 



' d_ 
dt 



v ± c (9" 
1 ± v c dz 



K± = 0, 



(55) 



where the so-called Riemann invariants are 



eo 



+ p(e') 
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Obviously, the Riemann invariants are constant, 
dK± 

c± 



dt 



0, 



(57) 



along the so-called characteristic curves C±(x,t), the positions x±(t) of which 
are determined by integrating 

dx± _ v ± c 

~df ~ r±^ • (58) 



The characteristics are the world-lines of sonic disturbances on the hydrody- 
namic flow pattern. For simple waves moving to the right (i.e., v > 0), one can 
prove [83] that 1Z + = const, and it suffices to consider the C_-characteristic. 
For simple waves moving to the left (v < 0), TZ- = const, and only the 
(^-characteristics have to be considered. 

For the initial conditions (53, 54) we determine the value of 1Z± at t — t , z = 
R± to be 1Z± = and thus 

eo 



If we consider the QGP EoS, p(e) 



c e 



Cq = const, this gives 



v(e)± = tanh y(e) p m 



tanh 



Vo =F 



Co 



1 + c^ 



In 



e 



tanht/o =F 



l_ (e / eo ) 2c o/(i+c§) 
l +(e/eo ) 2c o/(i+^) 



1 =F tanh y 



l-(e/e ) 2c ^ 1+c f 

l +(e/e0 )2 C 0/(l+cg) 



(60) 



This equation nicely shows the defining property of a simple wave, namely 
that there is a unique relationship between the value of the fluid velocity and 
its thermodynamic state [83]. Furthermore, with the help of (58) for C T , we 
can now calculate for any given e, < e < eo, the position z(t; e) at which 
this value for the energy density occurs for given t, 

Eq. (61) has similarity form, i.e., the profile of the rarefaction wave does 
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not change with time when plotted as a function of the similarity variable 

C = z t-R± 
■»± — t-t • 

Let us now complete the solution of the hydrodynamic problem. Causality 
requires that the initial conditions (53, 54) remain unchanged for \(±\ > 1. 
Thus, we only have to determine the solution in the range — 1 < (± < 1, i.e. 
in the forward light cone (remember that ( + (C_), obtained from the C_ (C + )- 
characteristic, eq. (61), corresponds to y + (y_), obtained from 1Z+ = const 
(1Z- = const) condition). 

For the forward going rarefaction wave, v > 0, generated on the right end of 
the initial streak, from eqs. (60, 61) we infer that the head of this wave (the 
point where the rarefaction of matter starts, i.e., where the energy density 
e starts to fall below e ) travels with the velocity tanh w-co to the left 0n 

u/ J 1— co tanh yo 

the other hand, the base of the rarefaction wave (the point where the vacuum 
ends, i.e., where e starts to acquire non- vanishing values) travels with light 
velocity v — 1 to the right. Thus, the energy density as a function of ( can be 
written as 



e(C+) = e • < 



1 



1 < /" < tanh yo— cp 

— ^+ — 1-co tanhyo 

1 + tanh yo 1 — Co 1 — C4 



1 — tanh y 1 + c 1 + C-4 

tanh yp— cp 
1— Co tanh yo 



tanhyp— cp < (^ + < X 



(l+ C 2)/2c 



(62) 



The velocity can then be inferred from (60), or simply from (61): 
tanhy(C+) = < 



tanh^ , -1<C + <SS§S 



. 1 + C+C ' i-^'t^W 



tanh yo— cq 



< C+ < 1 • 



(63) 



The similar analyzes for the backward going rarefaction wave, v < 0, generated 
on the left end of the initial streak, gives: 



e(C-) = e • < 



1 , 



tanhyo+cp a ^ j 
L+co tanhyo — ^~ — 



1+co tanhyo 

1 - tanh y 1 - c 1 + C- 
1 + tanh y 1 + c o 1 — C- 

_ y < /" < tanhyp+cp 

— 1+co tanhyo ' 



(l+cg)/2co 



(64) 
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tanhy(C-) 



tanh y 
C- ~ c 
I 1 - C-Co 



tanh yo+cp ^ /■ ^ ^ 
l+c () tanhj/o — S- — 

_ ^ ^ /■ < tanh yo+cp 
— l+c tanh?/o 



(65) 



This simple analytical solution is valid as long as the two rarefaction waves 
did not overlap in the middle of the system. Further evolution becomes more 
complicated and doesn't have a similarity form any more. 

The point where two rarefaction waves start to overlap, 
found from the system 

A / z j-\ — z-R+ _ tanh y -c Q 
=>+V > / t— to 1— eotanhj/o 



(66) 

t—<0 1+cotanhyo ' 



£ ^ z ^ _ z-R- _ tanhj/o+cp 



which gives 

R+(l - tanh y c ) (tanh y + c ) 
z<wer = — R_(l + tanh yoc ) (tanh y — Co) > (67) 
2c (l - tanh 2 y ) 

, _. . (zover - R+)(l - tanhy c ) , RQ , 
(tanh yo - c ) 



Thus the analytical solution (62-65) valid for t < t < t over . 

For the symmetric initial state, R + = — R- = R, eqs. (67,68) became 

tanhyo(l-cg) 
Co(l - tanh y ) 

1 — tanhy c o 

^troer = ^0 + 77 7 r~2 7 ■ (70) 

Cq(1 - tanh y ) 



Thus, based on the solution presented above we can have a more advanced 
description of the final streaks than the one presented in section 3.3. Let us 
assume that the homogenous final streaks, with some e/, yj, are formed in the 
CRF when the larger of initial streaks reach the rapidity yf. Up to this point 
the fluid cell trajectories are the same as those discussed in section 3.3, but 
after the homogeneous final streak is formed, it starts to expand according to 
the simple rarefaction wave solutions presented above. Figure 7 (B) shows new 
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trajectories of the streak ends, and Figures 7 (C,D) present energy density and 
rapidity profile for this expanding streak. 

Such an initial state with expanding streaks will also help to avoid the problem 
which may cause the development of numerical artifacts, namely a step-like in 
beam direction initial energy density distribution (output of ESRM): it has a 
jump from e(xj,j/j ) = const inside the matter to in the outside vacuum (of 
course, where is also a jump of e as a function of x, y, but it is much smoother 
and this is not the direction of the initial expansion). In order to avoid (or 
at least to suppress) the effect of this, it was proposed in Ref. [84] to smooth 
over initial energy density distribution, for example by a Gaussian shape. Our 
simple analytic solution smooths over this jump in a natural way. 

3. 6 Initial conditions for hydrodynamical calculations 

8r 
7- 



6- 




Fig. 10. Final rapidity, yj, profiles of the final streaks in the reaction plane for 
Au+Au collision at eo = 65 GeV/nucl, A = 0.08. The rapidities of the final streaks 
in CRF are calculated according to eq. (41). Our profiles are in agreement with the 
schematic sketch in paper [76] and with the scenario discussed in "firestreak" model 
(see Fig. 4). 

In this section we present the results of our calculations. We are interested 
in the shape of the QGP formed when the expansion of the combined target- 
projectile streaks stop and their matter is locally equilibrated. This will be the 
initial state for further hydrodynamical calculations. The time, th, at which 
we start the hydrodynamical description is a second (in addition to A) free 



34 



parameter of our model. Of course, th should be larger than the time of final 
streak formation, at least in the central, most hot and dense region. For the 
peripheral streaks the string tension is low, and the transparency is large, 
but peripheral matter does not play a leading role in further hydrodynamic 
expansion. Therefore, we will also build the final streaks (yf, ef) for peripheral 
streak-streak collisions, with lengths, Alf, corresponding to the lengths of the 
interacting region at the moment t = th, even if the final rapidity, yf, was 
not yet achieved for this particular collision.^ On the other hand the time, 
th, cannot be longer than the time for which the simple analytical solution 
presented in the previous section is valid. I.e., t h < Minimum{t f or m,i + tover,i} , 
where i/ rm,« is a time of formation of the final streak for particular streak- 
streak collisions i: when the larger of two initial streak gets rapidity yf (given 
by eq. (43) in Ref. [70]); t over>i is given by eq. (68). In this section we present 
all the results at th = 5 fm/c, while for e = 100 GeV/nucl, b = 0.25 • 2 Ra u , 
A = 0.08 the maximal possible th = Minimum{tf orm ^ +t over)i } ps 5.04 fm/c. 

The reader may notice that all of the above discussion about final streaks is for 
the energy density distribution. What about the baryon charge distribution in 
the final streaks? Our knowledge here is very poor and, in principle, we are free 
to assume anything between two extreme cases: A) the complete transparency 
picture, where charges would remain at the ends of the final streak; this is 
the qualitative picture in Ref. [71]; B) homogeneous in baryon charge final 
streaks, the simplified situation discussed in Refs. [68-70] 

First, we present the initial rapidities for final streaks for different x coor- 
dinates in the reaction plane - Fig. 10. Fig. 11 shows the energy density 
distribution in the reaction plane with non-expanding final streaks (Au+Au 
collision at e = 65 GeV/nucl, b = 0.25 • 2 R Au , A = 0.08). Fig. 12 shows the 
energy density distribution in the reaction plane for the same collision, but 
with expanding streaks. The QGP forms a tilted disk for b ^ 0. Thus, the 
direction of fastest expansion, the same as the largest pressure gradient, will 
be in the reaction plane, but will deviate from both the beam axis and the 
usual transverse flow direction. So, the new flow component, called "antiflow" 
or "third flow component" (see section 2.5), appears in addition to the usual 
transverse flow component in the reaction plane (for non-expanding streaks 
this was shown in [78-80]). With increasing beam energy the usual transverse 
flow is getting weaker, while this new flow component is strengthened. 

4 The calculations of the initial state presented in this section are done for col- 
liding nuclei with a homogeneous density distribution. The modification is pretty 
straightforward - we can use the Wood-Saxon density distribution for the Au nu- 
cleus in our case. This will lead to the increasing of the density in the middle of the 
colliding nuclei, therefore the middle steaks will contain more partons, will generate 
a stronger field and, thus, will form the final streaks faster. While the peripheral 
streaks will contain less matter, therefore their dynamics, which is oversimplified, 
will have less influence on the final state. 
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Fig. 11. Au+Au collision at e = 65 GeV/nud, b = 0.25 • 2 R Au , A = 0.08 (the 
parameter A was introduced in eq. (32)), E = T 00 is presented in the reaction plane 
as a function of x and z for different times in the laboratory frame. The final shape 
of the QGP volume is a tilted disk ~ 45°, and the direction of the fastest expansion 
will deviate from both the beam axis and the usual transverse flow direction and 
will generate the third flow component [78-80]. Note that the initial state for hydro 
looks pretty much like the one discussed in the "firestreak" model (see Fig. 4) , since 
the final streaks are moving with rapidities presented in Fig. 10. These calculations 
are done for the symmetrized T^ v (see eq. 35), but the results do not differ too 
much from what is presented in Ref. [70]. 

More results for full RHIC energy (e = 100 GeV/nucl) for different string 
tensions and impact parameters are presented in Figs. 13-14. 

Note that the initial state with non-expanding streaks (Figs. 11) looks pretty 
much like the one discussed in "firestreak" model (see Fig. 4), since the final 
streaks are moving with rapidities presented in Fig. 10. If we let final streaks 
expand, this smooths over picture, but most of the matter, nevertheless, keeps 
similar energy density profile and velocity distributions. 

From Figs. 11 and 13-14 we may see that for the central collisions at Eq = 
100 GeV/nucl the maximum energy density E = T 00 reaches as high as 
E ma x ~ 100 GeV/fm 3 . It seems to be larger than what one should ex- 
pect from the Bjorken estimate (see section 2.3) - eo = 4.3 GeV/fm 3 for 
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^/s = 130 A GeV [75]. This is natural since our initial energy density distri- 
bution is much more peaked as a function of z than is the one obtained in the 
Bjorken picture (our energy density distribution is also peaked as a function 
of transverse coordinate {x,y), but here the difference seems to be not that 
large, see for example [33]). One should, nevertheless, keep in mind that we 
have also a huge initial pressure, P = c^e, which will cause a fast hydrody- 
namical expansion, and thus, as discussed in section 2.3, will reduce ■ 

This is a usual situation for hydrodynamical calculation, for example, in Ref. 
[33] eo = 23 GeV j fm 3 has been used to reproduce the elliptic flow at RHIC for 
e = 65 GeV/nucl, although the Bjorken initial geometry had been assumed. 

In a recent work [77] another scenario of the creation of the third flow com- 
ponent based on the shadowing picture was proposed. The calculations were 
performed in the framework of QGSM. The shadowing effect manifests itself 
in the following way - hadrons emitted at small rapidities in the antiflow direc- 
tion can propagate freely, while hadrons emitted in the normal directed flow 
direction still remain within the expanding subsystem of interacting particles. 
We would like to point out how one can distinguish this scenario from the 
one discussed above in this paper. The tilted initial state has the maximal 
effect for non-zero, but not very large, impact parameters, approximately up 
to b = 0.5(^,4 + Rb), i-e., in semicentral reactions. For semiperipheral and 
peripheral collisions our string rope tension, a, would be small, and we would 
have to wait for a long time for the final streak formation, i.e., at the initial 
stages we would have one-dimensional Bjorken transparency. On the other 
hand, shadowing just becomes stronger with increasing impact parameter and 
gives the most drastic effect for semiperipheral and peripheral collisions. 



4 Conclusions 



A realistic and detailed description of an energetic heavy ion reaction re- 
quires a Multi Module Model, wherein the different stages of the reaction are 
each described with a suitable theoretical approach. The initial stages are the 
most problematic at RHIC energies. We tried to show that for the qualitative 
understanding of the basic features of the initial stages a phenomenological 
"firestreak" model can be used. For the detailed calculations of the initial state 
the Effective String Rope model has been developed. This model is based on 
earlier Coherent Yang-Mills field theoretical models and introduces an effec- 
tive string tension based on Monte-Carlo string cascade and parton cascade 
model results. Contrary to earlier expectations based on standard string ten- 
sions of 1 GeV/fm, which lead to the Bjorken model type of initial state, the 
effective string tension introduced in our model for collisions of energetic mas- 
sive heavy ions causes limited transparency. The increased string tension is a 
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consequence of collective effects related to QGP formation. These collective 
effects in central and semicentral collisions lead to an effective string tension of 
the order of 10 GeV/fm and consequently cause much less transparency than 
earlier estimates. The resulting initial locally equilibrated state of matter in 
semicentral collisions takes a rather unusual form (see Figs. 11, 12, 13-14) and 
the resulting flow velocity field (see Fig. 10, 7 (D)) is in agreement with the 
"firestreak" phenomenological picture. Such an initial state will manifest itself 
by the asymmetry of the produced collective flow, even after the subsequent 
hydrodynamical expansion [78-80]. 
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A Initial conditions after string creation 

Our conserved quantities are (36, 37) 

Q° = J f 00 dV = Ax Ay J f 00 dz , (A.l) 

Q 3 = J f 03 dV = Ax Ay J f 03 dz , (A.2) 



where T 00 and T 03 are given by eq. (38). Before string creation the initial 
values of the modified energy-momentum tensor, T^, are 

f ™ = f 2 °° = e cosh 2 y = 2 n m , (A.3) 



2 



T 2 Ud = -TT = e tanhy cosh 2 y = — n mv (A.4) 



where m is the nucleon mass, Eq is the initial energy per nucleon, and we have 

used cosh 2 = 7o = (m) ■ v ° = tanht/ is the initial velocity, vq = 1 is a 
good approximation for ultra-relativistic heavy ion collisions. Thus, 

Q° = Ax A y f— V n m(Zi + l 2 ) , (A.5) 
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2 



Q 6 = Ax A y ( — ) n m(Z 2 - Zi)u , (A.6) 



where l\ and Z 2 are the initial lengths of streaks (see Fig. 5), Ax, Ay are the 
grid sizes in x and y directions. Thus the rapidity of CM frame is 

tanh » = S = s- (A7) 



where 

f = 

^2 — ^1 



M = £±ii . (A.8) 



After string creation - 

T 00 = ei cosh 2 yi + Cgei sinh 2 y x + e 2 cosh 2 y 2 + c 2 ^ sinh 2 t/ 2 

+ ia 2 + B - ^- no e^e 2 ^ - ^n e~ yc 'e^ , (A.9) 



T 03 = ei(l + eg) cosht/x sinh?/! + e 2 (l + c 2 ,) coshy 2 sinh j/ 2 
^t^e^e 2 ^ + ^-n e~ yo e- 2y2 , (A. 10) 



At the point of complete penetration of streaks, t — t — (l 1 + l 2 )/2 (see Fig 5), 
we introduced energy densities ei(to) and e 2 (t )- We assumed transparency, 
i.e., that complete penetration happened so fast, that the field itself that was 
created during this time did not have time to stop partons. So, the rapidi- 
ties 2/1(2) (^o) — ~ l/o(l/o) correspondingly, and the proper baryon densities did 
not change, and thus, the baryon current is conserved. Then the energy and 
momentum conservation laws can be written in the form: 



Q° 



Ax Ay 



l + c )cosh y -c (ei(io)ii + e 2 (t )l 2 ) 



+ (^ + 5)(/ 1 + / 2 ) + ^f^(/ 2 + / 2 ) , (A.ll) 



Q 3 

Ax Ay 



;i + eg) 



cosh 2 



yo 



-ei(t )k + e 2 (t )k) 



an e 



-yo 



(l\ - f 2 ) 



(A.12) 
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We neglect C 2 , next to (1 + c 2 ,) cosh 2 y in eq. (A. 11), then eqs. (A. 11, A. 12) 
may be solved 

2 

/ , n n m \ + B h+l 2 an e- yo 

ei to =7—^-^2 77, 7 — T2 *i > A - 13 

1+Co fe) (l+c§) 2/1 4(a) (1+cl) 

2 

^°)-s- 7 V g ffi- rf yn fa- < a - 14 ) 
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Fig. 12. Simulation of the same collisions as in Fig. 11 (Au+Au, £o = 65 GeV/nucl, 
b = 0.25 • 2 Rau, A = 0.08), but with expanding final streaks. The middle subplot 
presents e(x, z) in the reaction plane for th = 5 fm/c. The upper subplot shows the 
e profiles for different x = const. We can clearly see three regions - two of forward 
and backward rarefaction and middle where the initial energy density, ej (eq. (42)), 
is still preserved. The lower subplot shows the e profiles for different z = const. 
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Fig. 13. Au+Au collision at e = 100 GeV/nucl, A = 0.065 , E = T 00 is pre- 
sented in the reaction plane as a function of x and z for = 5 fm/c. Subplot A) 
(b = 0.5 • 2 Rau)-, subplot B) (b = 0.25 • 2 Rau)- The QGP volume has a shape of a 
tilted disk and may produce a third flow component. 
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